# R-package of Sequential T-test Analysis of Regime Shifts (rSTARS)
# ref to: doi:10.1093/icesjms/fsz148
# ref to: https://github.com/LStirnimann/rstars
rstars <- function(data.timeseries = PDO, l.cutoff, 
                   pValue = 0.05, 
                   Huber = 1, 
                   SubsampleSize = (l + 1) / 3) {
  
  TS <- data.timeseries
  l <- l.cutoff
  Nsub <- SubsampleSize
  
  
  DT = 0
  RSI_mat = matrix(0, nrow = length(TS[, 1]), length(TS[1, ]))
  RMean_mat = matrix(0, nrow = length(TS[, 1]), length(TS[1, ]))
  
  
  #### attaching the data set and removing of red noise ####
  for (TIMESERIESindex in 2:length(TS[1, ]))
  {
    
    X = ts(TS[, TIMESERIESindex])
    N = length(X)
    
    if (N < l)
    {
      stop("CutOff cannot be > Time series length")
    }
    
    if (Nsub > N)
    {
      Nsub = N
    }
    
    
    #===================#
    ####  STARS 3.2  ####
    #===================#
    
    #freedom degree
    df = 2 * l - 2
    
    #two tailed test
    t_stu = abs(qt(pValue / 2, df))
    
    #Variance and Sigma calcualation for DIFF formula
    A = var(X[1:l])
    
    for (i in 2:(length(X) - l + 1))
    {
      B = var(X[i:(i + l - 1)])
      A = rbind(A, B)
    }
    
    #Sigma square
    Sigma_s = mean(A)
    
    #between mean values of two subsequent regimes that would be statistically
    #significant according to the Student’s t-test
    diff = t_stu * sqrt((2 * Sigma_s) / l)
    
    #====================#
    #     core steps     #
    #====================#
    
    
    vRMean = 0
    RSI = seq(0, 0, length.out = length(X))
    
    R1 = X[1:l]
    RegimeMean = WeightedAverage(R1,Sigma_s,Huber)
    changepoint = 1
    n1 = 0
    
    for (intYear in 2:length(X))
    {
      
      if (is.na(RegimeMean) || RegimeMean == '')
      {
        break
      }
      
      if (X[intYear] > (RegimeMean + diff))
      {
        sumofWeights = 0
        cusumUP = 0
        Xdev = 0
        for (t in intYear:(intYear + l - 1))
        {
          if (t > length(X))
          {
            if (sumofWeights > 0)
            {
              break
            }
          }
          
          Xdev = (X[t] - RegimeMean - diff) / sqrt(Sigma_s)
          
          #determine the weight of the normalized deviation
          if (Xdev == 0)
          {
            Xweight = 1
          }
          
          else if (Xdev != 0)
          {
            Xweight = min(1, (Huber / abs(Xdev)))
          }
          
          #sum weights and weighed values
          sumofWeights = sumofWeights + Xweight
          cusumUP = cusumUP + (Xdev * Xweight)
          
          #check if cusum turns zero
          if (cusumUP < 0)
          {
            cusumUP = 0
            break
          }
        }
        cusumUP = cusumUP / sumofWeights
        
        RSI[intYear] = cusumUP
      }
      
      else if (X[intYear] < (RegimeMean - diff))
      {
        sumofWeights = 0
        cusumDown = 0
        Xdev = 0
        for (t in intYear:(intYear + l - 1))
        {
          if (t > length(X))
          {
            if (sumofWeights > 0)
            {
              break
            }
          }
          
          Xdev = (X[t] - RegimeMean + diff) / sqrt(Sigma_s)
          #determine the weight of the normalized deviation
          if (Xdev == 0)
          {
            Xweight = 1
          }
          else if (Xdev != 0)
          {
            Xweight = min(1, (Huber / abs(Xdev)))
          }
          
          #sum weights and weighed values
          sumofWeights = sumofWeights + Xweight
          cusumDown = cusumDown + (Xdev * Xweight)
          
          #check if cusum turns zero
          if (cusumDown > 0)
          {
            cusumDown = 0
            break
          }
        }
        cusumDown = cusumDown / sumofWeights
        RSI[intYear] = cusumDown
      }
      
      
      else if (RegimeMean - diff <= X[intYear] & X[intYear] <= RegimeMean + diff)
      {
        RSI[intYear] = 0
      }
      
      #check for the situation when the test is not over for the last
      #change point, but we are too close to the end of the time series
      if (abs(RSI[intYear] > 0 & intYear > (length(X) - l + 1)))
      {
        break
      }
      #------------------------------------------------------------------#
      
      if (RSI[intYear] == 0)
        #intYear is not a new changepoint
      {
        if ((changepoint + l) <= intYear)
        {
          #recalculate regime mean and Diff
          #currently Diff remains constant for the entire process /series
          n1 = intYear - changepoint + 1
          for (n in 1:n1)
          {
            R1[n] = X[changepoint + n - 1]
          }
          RegimeMean = WeightedAverage(R1,Sigma_s,Huber)
        }
      }
      
      
      if (RSI[intYear] != 0)
        #regime shift is detected
        #intYear is a new changepoint
      {
        changepoint = intYear
        #recalculate regime mean and Diff
        #currently Diff remains constant for the entire process /series}
        R1 = 0
        for (n in 1:l)
        {
          R1[n] = X[changepoint + n - 1]
        }
        RegimeMean = WeightedAverage(R1,Sigma_s,Huber)
        
      }
    }
    
    X1 = TS[, TIMESERIESindex]
    S = 1
    for (i in 1:length(RSI))
    {
      if (RSI[i] != 0)
      {
        E = (i - 1)
        MeanRegime = WeightedAverage(X1[S:E],Sigma_s,Huber)
        vRMean1 = rep(MeanRegime, length(X1[S:E]))
        vRMean = c(vRMean, vRMean1)
        S = i
      }
      if (i == length(RSI))
      {
        E = (length(RSI))
        MeanRegime = WeightedAverage(X1[S:E],Sigma_s,Huber)
        vRMean1 = rep(MeanRegime, length(X1[S:E]))
        vRMean = c(vRMean, vRMean1)
      }
    }
    
    vRMean = vRMean[-1]
    RSI_mat[, TIMESERIESindex] = RSI
    RMean_mat[, TIMESERIESindex] = vRMean
    
  }
  
  ####Saving tables of regimes avarege (tsMean.txt), RSI.txt and Filtered time series (Filredts.txt)####
  
  colnames(RMean_mat) <- colnames(data.timeseries)
  colnames(RSI_mat) <- colnames(data.timeseries)
  
  RMean_mat[,1] <- data.timeseries[,1]
  RSI_mat[,1] <- data.timeseries[,1]
  
  return(list(
    mean = RMean_mat,
    rsi = RSI_mat
  ))
}

#function for Average Mean
#EstAve = estimated average of the range
#Xdev = deviation from the median
#Xweight = weights

WeightedAverage <- function(W,Sigma_s,Huber)
{
  EstAve = mean(W)
  for (k in 1:2)
  {
    SumofWeights = 0
    SumAve = 0
    for (i in 1:length(W))
    {
      Xdev = ((W[i] - EstAve) / sqrt(Sigma_s))
      #determine the weight of the normalized deviation
      if (is.na(Xdev) || Xdev == '')
      {
        break
      }
      if (Xdev == 0)
      {
        Xweight = 1
      }
      else if (Xdev != 0)
      {
        Xweight = min(1,(Huber / abs(Xdev)))
      }
      
      
      #sum weights and weighed values
      SumofWeights = SumofWeights + Xweight
      SumAve = SumAve + Xdev * Xweight
    }
    SumAve = SumAve / SumofWeights
    SumAve = SumAve * sqrt(Sigma_s) + EstAve
    EstAve = SumAve
  }
  WeightedAverage = EstAve
}
